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Введение. Целью работы является разработка нового универ- 
сального численного метода решения краевых задач для ли- 
нейных уравнений эллиптического типа. 

Материалы и методы. Предложенный метод основан на 
приведении исходного уравнения математической физики к 
более простому неоднородному уравнению с известным фун- 
даментальным решением. Затем предлагается произвести 
переход к неоднородному интегральному уравнению с ядром, 
выражаемым через известное фундаментальное решение. По- 
лученное интегральное уравнение совместно с граничными 
условиями решается численно. В результате было получено 
искомое приближенное решение, потенциал поля в аналити- 
ческом виде, что позволило не только находить приближен- 
ное значение потенциала поля в любой точке области реше- 
ния, но и дифференцировать этот потенциал, причем без за- 
метной потери точности. Это свойство разработанного чис- 
ленного метода выгодно отличает его от традиционных чис- 
ленных методов решения краевых задач, таких, например, как 
метод конечных элементов. 

Результаты исследования. Для подтверждения эффективно- 
сти предложенного численного метода решена двумерная и 
трехмерная краевые задачи с известными решениями. Полу- 
чены зависимости погрешности численного решения от числа 
линейных уравнений в результирующей системе. Показано, 
что даже при небольшом числе уравнений в системе (порядка 
нескольких сотен) достигается точность решения на уровне 
сотых долей процента. Еще одной важной иллюстрацией эф- 
фективности предложенного метода является решение кван- 
тово-механической задачи для одномерного и двумерного 
квантового осциллятора. Показано, что рассматриваемый 
метод позволяет находить собственные значения энергии и 
собственные функции с приемлемой точностью. Разработан- 
ный численный метод позволяет существенно расширить 
область применения традиционного метода точечных источ- 
ников поля при решении прикладных задач моделирования 
полей различной физической природы, включая задачи на 
собственные значения. 

Обсуждение и заключения. Полученные результаты подтвер- 
ждают, что физическое поле, описываемое практически лю- 
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бым линейным уравнением эллиптического типа, можно 
представить в виде суперпозиции полей точечных источни- 
ков, удовлетворяющих более простому уравнению, решение 
которого находится с помощью метода точечных источников 
поля. Таким образом, представленный в данной работе чис- 
ленный метод можно рассматривать как обобщенный метод 
точечных источников поля. 


Ключевые слова: фундаментальное решение, метод фунда- 
ментальных решений, метод точечных источников, уравнения 
эллиптического типа, краевая задача. 
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Введение. Основными методами численного решения краевых задач для уравнений математической физики в 
настоящее время являются метод конечных элементов (МКЭ) и, в значительно меньшей степени, метод конечных раз- 
ностей (МКР). Эти методы могут использоваться при решении практически любых линейных уравнений, что позволя- 
ет считать МКЭ и МКР универсальными методами решения краевых задач. Область применения других методов 
обычно ограничивается вполне определенными типами уравнений математической физики. Тем не менее, методам 
МКЭ и МКР присущи определенные недостатки. Это, например, недостаточно быстрая сходимость численного реше- 
ния с уменьшением шага разностной сетки и обусловленная этим высокая погрешность результата. Определенные 
трудности возникают, если граничные условия содержат производную по нормали от искомой функции. Если кроме 
искомой функции необходимо определить и ее производные, то для этого необходимо производить численное диффе- 
ренцирование, что также ведет к появлению дополнительной, весьма значительной погрешности. В связи с этим про- 
блема поиска новых численных методов решения краевых задач продолжает оставаться актуальной. 

Одним из эффективных методов решения краевых задач, который может являться альтернативой МКЭ и 
МКР, является метод точечных источников поля (МТИ) [1-5]. Этот метод может использоваться для решения широ- 
кого круга задач математической физики. Наиболее эффективно использование этого метода при решении краевых 
задач для уравнений эллиптического типа: Лапласа, Гельмгольца, бигармонического уравнения [6-12]. Имеются све- 
дения о возможности и эффективности использования этого метода при решении краевых задач для уравнений пара- 
болического типа, для волновых уравнений [13-15], при моделировании полей упругих напряжений в деформирован- 
ных твердых телах [16-20] ит. д. 

Для МТИ характерны высокое быстродействие в сочетании с малой вычислительной погрешностью, а также 
простота компьютерной реализации. Для возможности использования МТИ необходимо располагать фундаменталь- 
ными решениями исходного уравнения в аналитическом виде. Однако при решении некоторых задач математической 
физики фундаментальные решения неизвестны, что является препятствием для использования МТИ. Тем не менее, 
как показано в [21-22], имеется возможность получения фундаментальных решений численно, с последующим их 
применением в МТИ. С помощью подхода, используемого в [23-24], искомое решение краевой задачи было получено 
без предварительных вычислений фундаментальных решений. Это позволило значительно расширить область приме- 
нения МТИ. 


Интегральные уравнения для численного решения краевой задачи. Пусть в области О фундаментальное 


решение $ (х, К) уравнения [И =0 известно. Требуется найти решение краевой задачи для неоднородного уравнения 


(Е+К(к))И (+) = Г (+), гео, (1) 
с условиями на границе 00: 
(ги (г) =$(‚), гебО, (2) 
где Г, Кг) , А(г) — линейные операторы в области ®© и на границе 00. 
Для определенности будем считать, что операторы Ц») и А(г) имеют вид: 


(г) =а(г)+В(г)-У ‚ ГЕО, 
(г) (+8). ГЕДО. 


Здесь а(к), В (+) } @(х), В(х) — заданные скалярные и векторная функции в области © и на границе 00, 


\У — оператор Гамильтона. 
Ищем решение задачи (1)-(2) в виде: 


0 (›) =и(г)+ъ(к) р 
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где функция у (г) удовлетворяет однородному уравнению 
м=0, (3) 
а неизвестная функции и (^) есть частное решение уравнения 
Ги = } —Цг)и —Цг)у(г) р 
правую часть которого }-—/ (г)и —1 (*) У= р(*) можно рассматривать как плотность р(») некоторого фиктивного 
заряда, распределенного в области о. Отсюда ДЛЯ функции и (г) получаем интегральное уравнение 
и(")=[8 ("кур (К)ао, 
ИЛИ 
и(")= [8 (к) (л(®)-Цк)у(к)-Кк)и(к)) ао. (4) 


Для функции у(г) ‚ как в МТИ, используется приближенное выражение 


о (= ив (ри) (5) 


Здесь №, — число точечных зарядов 4», моделирующих поле, располагаемых в точках с координатами р» за пределами 
области решения задачи ©, вблизи ее границы. 
Теперь уравнение (4) принимает вид: 


ие) ав (®-К®) (вр) Ви) 40 © 


Если оператор КК) представляет собой скалярную функцию (это справедливо, если решается, например, 
стационарное уравнение Шредингера), то, вычисляя интеграл в правой части соотношения (6) численно, можно, как 


будет показано ниже, построить алгоритм для нахождения значений функции и(г) в узловых точках области ©. В 
более общем случае подействуем на уравнение (6) оператором [(г). Введя при этом новую функцию (г) = [(г)и(к), 


получим для этой неизвестной функции интегральное уравнение: 
ве) ав (вн) -К®) Уч (Вл, ) Чо. (п) 
Решая уравнение (7) численно находим функцию »’(г) ‚ после чего, используя уравнение (6), записанное в виде 
ие) = [85 ®)[ (®- ЦК) ав ( Вр) (8) 46, ®) 
находим искомую функцию и(›). 


Пусть теперь г Е ОО. Подействуем на уравнение (8) при гЕДО оператором А(г) | 


А(г)и(к) = Ав (вв (В) - КЮ) ( В) 40». 
Теперь граничное условие (2) можно записать в виде: 
Абв (вв (В) - КУ (Вл, )-"(®)) ао, ал) ) = ф(+), ГЕДО. (9) 


Соотношения (7)-(9) можно использовать для численного решения исходной краевой задачи. Построим алгоритм ре- 
шения этой задачи. 


Алгоритм численного решения краевой задачи. Область решения © разбивается на элементарные участки 
площадью (объемом), равным с2;. Число таких участков обозначим как №,. Внутри каждого участка выбираем точку, 
узел с координатой К;. Некоторые из точек, узлов г„=А;, число которых М№с=М№., располагаются на границе области 
решения. Теперь уравнение (7) принимает вид: 


ве) ов) (в), = > [ (в) Ха) (Вр) Ков ( в), (10) 


Это уравнение должно выполняться во всех №, узловых точках, включая М№с граничных узлов. Запишем уравнение 
(10) для каждого узла К;. Получим: 


(в) + (вв (В.В, ( в), = ХЦВ)в(в,в, [ 78, Ха(® (Вл, и) 
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В результате получаем М№, уравнений для М№М,-+№ с неизвестных. Это №, значений функции (г) в №, узловых 


точках №, и №. =М№ с зарядов 4 к, располагаемых в точках с координатами р». Для получения замкнутой системы урав- 
нений необходимо дополнить систему (11) уравнениями, учитывающими граничные условия для решаемой задачи. 
Для этого используем интегральное уравнение (9). 

На границе ОО, в точках с координатами г„ располагаем граничные узловые точки, количество которых равно 
М№.=М№ с. Заменяем приближенно интегралы в (9) интегральными суммами, как это сделано при получении соотноше- 
ний (10), (11). В результате для каждого граничного узла получим уравнение: 


ХА )в (вы, 7(в,)- Хи (®,) (Вр) (+ Ада (ри) = 465) (12) 


Решая объединенную систему уравнений (11), (12) находим М, значений функции и(ь) в узловых точках об- 
ласти ©, и № с значений точечных зарядов 4х. Затем с помощью соотношения (8) вычисляем приближенное значение 


функции и(г) в произвольной точке области ©, и, обозначив "(К | как и, а Г (в,) как ],, получаем искомое ре- 


шение исходной задачи в виде соотношения: 


№ № № 
(г) = Х (вк, К, } 4.8 (К›в, =) м + и (т.р, ). (3) 
Следует обратить внимание на то, что соотношение (13) дает приближенное аналитическое выражение для 
искомой функции Ц (*) . Это значит, что с полученным решением И (^) можно поступать как с любым другим анали- 
тическим выражением. Его можно дифференцировать, интегрировать, производить другие действия, и при этом не 
возникает дополнительная численная погрешность. Например, частную производную по координате х от функции 
и (г) можно вычислить с помощью формулы: 
бИ (г) м. д8 (ьк,) м м бв(р,) 

=У ( ‚КЕ, ) ха, 8(к, —, + А (14) 

а. (Е) из (Вр, №, в) + 5.4, 5 


Аналогичным образом можно получить аналитические выражения для вычисления результатов других мате- 





матических операций над функцией И (^) : 


Решение тестовых задач. Ниже приводятся результаты решения тестовых задач, подтверждающих эффек- 
тивность рассматриваемого метода при решении краевых задач для уравнений эллиптического типа. Сначала реша- 
лась двумерная краевая задача. В области ©, представляющей собой квадрат со стороной [=2, центр которого совме- 
щен с началом координат, решалась задача для уравнения 


0 (") 90 (г) 





АО (г) +2" — бу + т (ху)0 (г) = У (г) ‚ ГЕО (15) 
с условиями на границе 0@ 
Е =$(+), гео. (16) 
4 бп 


Функции } (г) и (+) в соотношениях (15), (16) подбирались таким образом, чтобы они соответствовали 
точному решению краевой задачи (15)—(16) в виде: 

О (г) =0(х, у) =х' +у+зщ(ху). (17) 

При решении краевой задачи оператору Г. в уравнении (1) соответствовал двумерный оператор Лапласа Г=А 


и, соответственно, фундаментальное решение $ (г,К) имело вид: 
8(",®)= ЧР шр- |. 
2п 


р 


Оператору [(к) в (1) соответствовал оператор [(г) =2у =. +х? + зт (ху) ‚ а оператору Л (и) в (2) — оператор 


х у 
^(*)=1+ а 
4 бп 
При численном решении краевой задачи (15)-(16) задавалось число М, определяющее равномерный шаг й 
сетки в области решения (9. Точечные заряды 4%, моделирующие поле, равномерно располагались по периметру квад- 
рата со стороной [..=АЁ. Параметр А>1 определял удаленность зарядов от границ области решения. При этом слишком 
большое и слишком малое удаление зарядов, как и в стандартном варианте МТИ, нежелательны, так как это ведет к 
снижению точности решения задачи. Авторами при вычислениях использовалось К=1,3. Полное число узлов в области 
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2 

© равнялось №, =М№`, а количество точечных зарядов, моделирующих поле №.=4(М№-1). Таким образом, полное число 
>. 2 

уравнений в системе (11)—(12) составило №, =М№,+М№М.= № +4(М№-—1. 

После нахождения приближенного численного решения краевой задачи И», (г) ‚ при заданном полном числе 
уравнений Л =М, в системе (11)-(12), оценивалась максимальная относительная погрешность решения &вх. Для этого 
в области решения задачи случайным образом задавались координаты #1; точек, число которых полагалось равным 
М№.=100, после чего относительная погрешность вычислялась с помощью формулы 

тах |» (=)-0 (= ) 
8 тах = и 


тах 





> 


где ПО „„— максимальное по абсолютной величине значение точного решения (17) в области ©. Аналогичным обра- 


зом, наряду с максимальной относительной погрешностью оценивалась и среднеквадратичная погрешность. 
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Рис. 1. Зависимость максимальной относительной погрешности 
численного решения от количества узлов при решении двумерной и 
трехмерной краевой задачи 
Но. 1. Оерепдепсе о{ тахипит ге]аНуе еггог оЁ патепса] зоаНоп оЁ фе питбег 
оЁ подез ипаег зо]утз 6мо-птепз!опа! ап4 Фе #гее-Айпепз1опа! Боипдагу 
уаше рго ет 


На рис. | представлена зависимость максимальной относительной погрешности от полного числа уравнений 
№=М, в системе (11)-(12) (жирная сплошная линия). С ростом числа М от начального значения М=17 до значения 
М№=721 погрешность решения уменьшается более чем на два порядка, от &ик=2, 1107! до к=7,5-10—^. 

Для получения системы (11)-(12) использовался простейший метод численного интегрирования, аналог мето- 
да прямоугольников. Можно предположить, что использование более точного метода численного интегрирования 
позволит повысить точность решения краевой задачи. Действительно, использование метода Гаусса с третьим поряд- 
ком интегрирования позволило повысить точность численного решения почти на порядок. Это проиллюстрировано 
рис. 1, где тонкая сплошная линия представляет зависимость погрешности решения краевой задачи от числа уравне- 
ний М, рассчитанной при использовании метода Гаусса. 

Далее решалась трехмерная задача Дирихле. В области ©, представляющей собой куб со стороной [=2, центр 
которого совмещен с началом координат, решалась задача для уравнения 
ЗЕЕ). а 90 (г) 90 (г) 

+х + 
дх ду 0 





АИ (+25 +51 (>) (г) =/(*), гео. (18) 


Функции / (*) в уравнении (18) и в правой части условия Дирихле подбирались таким образом, чтобы они соответ- 


ствовали точному решению краевой задачи в виде: 
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О (г) =0 (х, у, =) =х' +у+зщт(ху)- =. (19) 
При решении краевой задачи оператору Г, в уравнении (1) соответствовал трехмерный оператор Лапласа [=А 


и, соответственно, фундаментальное решение & (® К) имело вид 





1 
Е 
Е 
20 ›0 0 . 
Оператору [(г) в (1) соответствовал оператор [(г)=2у - +х г ы- +91 (ху), а оператору Л(и) в (2) соответ- 
х У д 


ствовал оператор Л(г)=1. 


Пунктирная кривая на рис. 1 представляет зависимость относительной погрешности решаемой задачи Дирих- 
ле от размерности М, т. е. от числа уравнений в системе (11)-(12), для решения которой использовался аналог метода 
прямоугольников. Видно, что уже при относительно небольшой размерности численной задачи №=800 погрешность 
решения составляет =их=2,0-10. Можно предположить, что и в этом случае использование более точного метода 
численного интегрирования позволит снизить погрешность решения краевой задачи. 


Решение стационарного уравнения Шредингера. Фундаментальное значение при решении квантово- 
механических задач имеет уравнение Шредингера [25]. Аналитические решения этого уравнения могут быть получе- 
ны лишь для весьма ограниченного круга задач, преимущественно одномерных. Несмотря на широкий спектр имею- 
щихся численных методов решения уравнения Шредингера, проблема эффективных способов нахождения собствен- 
ных энергий и собственных функций для основного уравнения квантовой механики, особенно при решении много- 
мерных задач, продолжает оставаться актуальной [23-24]. Предложенный численный метод может, в ряде случаев, 
использоваться при решении квантово-механических задач. 


Пусть частица массой т совершает финитное движение в силовом поле У(!). Запишем стационарное урав- 

нение Шредингера в виде: 
дут (и). (20) 

Необходимо, решая уравнение (20), найти собственные значения энергии Ё и соответствующие им собственные функ- 
ции, волновые функции частицы (г) ; 

При решении краевой задачи для уравнения Шредингера оператору Ё в уравнении (1) соответствует оператор 
Лапласа [=Л, а оператор (г) в (1) есть скалярная функция /(г) = ТЕ - ту (к). 

При реализации алгоритма решения уравнения Шредингера сначала задается энергия частицы Е. Область ре- 
шения для уравнения Шредингера © берется близкой, но несколько большей классически доступной для частицы об- 
ласти движения, границы которой определяются уравнением 


У(^)=Е. 2 


Размеры области © подбираются таким образом, чтобы на ее границах волновая функция с заданной точно- 
стью равнялась нулю. Для получения нетривиального решения задается некоторая точка М внутри области ©, напри- 
мер, на границе классически допустимой области движения, в которой значение волновой функции заведомо не рав- 


няется нулю. В этой точке значение волновой функции полагается равным единице: (М) =1. После получения ре- 


шения волновую функцию можно легко перенормировать. 
После численного решения задачи описанным выше способом проверяется найденное значение волновой 


функции в точке С. Если окажется, что с достаточной точностью выполняется условие у(С) =0 › то это означает, что 
заданное значение энергии Е является собственным значением, а найденная волновая функция (г) является соб- 


ственной функцией для данной энергии ЕЁ. В противном случае подбираются такие два значения энергии Е! и Ё>> Е! 
при которых у(С) имеет различные знаки. После этого собственное значение энергии в интервале (ЁЕ1, Е) и соб- 
ственная функция находятся, например, с помощью метода половинного деления. Иллюстрацией такого подхода яв- 
ляется рис. 2, на котором представлен график зависимости функции $1п (%(с)) от значения квантового числа и для 


одномерного квантового осциллятора. 
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Рис. 2. Зависимость знака волновой функции в граничной точке С от значения 
квантового числа п 


Но. 2. Оерепдепсе оЁ \уауе ГапсНоп $121 аё Боип4дагу рошё С оп Фе уаше оЁ апапат питаЬег и 


Как известно, собственные значения одномерного квантового осциллятора определяются выражением 


1 
Е = пы › Где ® — частота осциллятора, а квантовое число п принимает целочисленные значения, начиная с 


нуля. Как продемонстрировано на рис. 2, вблизи целочисленных значений и знак у(С) изменяется с положительного 
(при меньших п) на отрицательное (при больших п). Вблизи такого перехода можно подобрать значение п=их, при 
котором \(С) с приемлемой точностью равняется нулю. Величина Ап отклонения пм от целочисленного значения 


определяет погрешность для численного значения собственной энергии. 


На рис. 2 также изображены переходы от отрицательных значений \(С) к положительным. По мере при- 


ближения к такому переходу значение \(С) резко возрастает по абсолютной величине, но не обращается в ноль. По- 


этому такие переходы не соответствуют собственным значениям энергии. 

Точность в определении собственных значений энергии увеличивается с ростом размерности М системы ли- 
нейных уравнений, к которой приводится численная задача для уравнения Шредингера. Наибольшая точность, есте- 
ственно, наблюдается при решении одномерной задачи Шредингера. Например, относительная погрешность ДЁ при 
вычислении первого собственного значения энергии одномерного осциллятора изменялась от значения 0,15 при №=20 
до значения 0,00025 при М=120. При решении двумерных и трехмерных задач точность вычислений значительно сни- 
жается. На рис. 3 представлены зависимости относительной погрешность ДЁ для собственных значений энергии дву- 
мерного осциллятора от размерности № численной задачи. Жирная сплошная линия на рис. 3 соответствует основно- 
му, нулевому энергетическому уровню, тонкая сплошная линия — второму энергетическому уровню, пунктирная ли- 
ния — четвертому энергетическому уровню. При решении двумерной задачи относительная погрешность порядка 
нескольких процентов легко достижима. Можно ожидать, что использование более точного метода численного инте- 
грирования, как и в предыдущем случае, позволит повысить точность численного решения квантово-механической 


задачи. 
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Рис. 3. Зависимость относительной погрешности от размерности системы 
уравнений для трех собственных значений энергии двумерного осциллятора 


Е. 3. Оерепдепсе оф ге]апуе еггог оп Читепз1опаШу оЁ зу$ет о{ едиаНопз 
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Описанный алгоритм численного решения стационарных квантово-механических задач не зависит от их раз- 
мерности. Решение может несколько усложниться при наличии вырождения энергетических уровней, когда одному 
собственному значению энергии соответствует несколько собственных функций. Например, при решении задачи для 
двумерного квантового осциллятора собственные значения энергии определяются двумя квантовыми числами п, ип,: 


Е=Е, +Е, = по (п, +п, + 1) . При численном решении краевой задачи область решения ©, в соответствии с (21), зада- 


ется в виде прямоугольника, несколько большего классически доступной для частицы области движения, с соотноше- 


х 


Г, о ы 
нием сторон —. = . Отсюда следует, что, например, первый энергетический уровень с п=п, +п, =1 два- 





жды вырожден. Этому энергетическому уровню соответствуют две собственные функции, рассчитываемые при 








Г 
п, =Ь М, =0 ип, =0, п, =1, и при соотношении сторон а = Зи 


У 


Г. 1 
р = = При другом соотношении сторон по- 


лучается смешанное состояние. 


Заключение. Полученные результаты показывают, что физическое поле, описываемое практически любым 


линейным уравнением эллиптического типа (1+1 (,)) 6 (г) = 1(г) можно представить в виде суперпозиции полей 


точечных источников, удовлетворяющих более простому уравнению (г) = Х (^) › решение которого находится с 


помощью МТИ. Поэтому разработанный численный метод можно рассматривать как обобщенный метод точечных 
источников поля (ОМТИ). Решение тестовых задач подтверждает эффективность ОМТИ при решении достаточно 
сложных краевых задач, включая трехмерные задачи и задачи на собственные значения, в том числе квантово- 
механические. Существенным преимуществом ОМТИ по сравнению с традиционными методами численного решения 
краевых задач, таких как МКР или МКЭ, является возможность получения решения в аналитическом виде, что позво- 
ляет производить с этим решением соответствующие математические операции без потери точности. Численные экс- 
перименты показали, что погрешность для производных по координатам имеет тот же порядок точности, что и по- 
грешность для искомого потенциала. 
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